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technique  is  applied  to  two  state  variable  constraint  problems,  in  one  of 
which  a penalty  function  is  employed  to  convert  the  constraint  problem 
to  an  unconstrained  one  in  addition  to  the  approach  considering  the 
constraints  directly.  For  this  same  problem  the  method  of  steepest 
descent  also  is  studied,  and  comparison  of  the  results  obtained  is  made 
and  discussed. 
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PREFACE 


The  Multiple  ere  trajectory  optimization  is  one  which  constantly 
confronts  Air  Force  flight  vehicle  systems,  both  in  air-to-air  and 
air- to  ground  operations.  The  most  difficult  of  problems  here  includes 
both  bounds  on  state  and  control,  and  yet  it  cannot  be  avoided  because 
this  is,  in  fact,  the  situation  in  Air  Force  flight  vechiles.  This 
report  appears  to  represent  one  of  the  most  important  pieces  of  work 
presenting  results  of  importance  both  for  flight  control  because  of  the 
greatly  effecient  algorithms  developed  in  this  report. 
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INTRODUCTION 
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The  past  decade  has  seen  considerable  progress  in  techniques  for 
optimisation  of  nonlinear  dynamical  systems.  The  development  of  large 
digital  computers  coupled  with  the  Interest  in  optimal  control  theory, 
particularly  in  optimising  of  spacecraft  trajectories,  has  inspired  s 
large  volume  of  literature  devoted  to  both  the  mathematical  theory  of 
optimal  processes  and  the  methods  for  obtaining  solutions  to  these 
problems.  Nevertheless,  from  the  computational  standpoint  the  class 
of  control  problems  with  constraint  state  variables  has  scarcely  been 
considered,  although  these  types  of  problems  often  occur  in  engineering 
practice.  Tor  example,  the  velocity  of  a vehicle  may  be  limited  by 
structure  breakdown  or  a motor  may  be  overloaded  to  prevent  safety  and 
reliability  of  operation.  Bryson,  Denham  and  Dreyfus  Cl],  [23*  and 
Starr  [33  have  treated  this  class  of  problems  using  the  steepest 
descent  technique  and  a suitable  combination  of  various  non-gradient 
techniques,  respectively.  Others  [23,  C**3  have  reduced  the  constraint 
problem  to  unconstrained  status  by  introducing  the  penalty  function  in 
place  of  the  constraints  on  the  state  variables. 

The  method  of  steepest  descent  is  excellent  for  finding  an 
approximate  solution  quickly,  but  it  often  exhibits  very  slow 
convergence,  whereas  other  techniques  frequently  face  the  problem  of 
computational  stability  in  the  solution  of  the  two-point  boundary  value 
problem.  It  is  hoped  that  the  method  of  conjugate  gradients  would 
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offer  an  improved  end  more  efficient  computational  aethod,  which  is  the 
objective  of  this  study. 

For  completeness,  some  basic  concepts  froa  functional  analysis  to 
be  used  in  the  sequel  are  given  in  Section  II,  and  a review  of  the 
conjugate  gradient  aethod  of  Hestenea  and  Stiefel  for  linear  operator 
equations  and  the  extension  to  nonlinear  operator  equations  are  given 
in  Section  III  with  emphasis  toward  control  applications.  Section  IV 
discusses  the  class  of  control  problems  to  be  considered.  The 
computational  aspect  of  the  state  variable  constraint  control  problems 
is  presented  in  Section  V.  The  application  of  the  method  of  conjugate 
gradients  to  this  class  of  optimal  control  problems  is  discussed.  An 
algorithm  is  given  showing  the  construction  of  the  sequence  of  control 
functions  that  extremize  a given  performance  functional.  Sections  VI 
and  VII  consider  two  practical  engineering  applications:  i)  a 

minimum  time  problem  in  two  dimensions  with  the  constraint  on  the 
state  variables  being  a circle,  and  2)  a re-entry  vehicle  problem 
with  altitude  constraint  for  which  the  range  is  to  be  maximized.  A 
comparison  of  the  rate  of  convergence  with  the  method  of  steepest 
descent  is  given  in  the  first  problem.  The  results  showed  that  the 
method  of  conjugate  gradients  provided  a higher  rate  of  convergence, 
but  not  as  rapid  as  for  the  cases  without  state  variable  constraint. 


SECTION  II 

BASIC  CONCEPTS  PROM  FUNCTIONAL  ANALYSIS 

Some  definitions  and  fundamental  theorems  from  analysis  which  will 
be  used  in  the  following  discussion  are  given  below.  Let  V be  a real 
Hilbert  space  with  inner  product  denoted  by  (•,*>.  u,  h are  elements 
in  V.  A function  t(h)  is  written  o(  | |h|  | ) if  — 0 as  | |hl  I 0, 

and  a function  t(h)  is  written  0(l|h||)  if  U(h)|  < N||h|l  as  llhl I ♦ 0 
where  N is  a positive  constant  and  llhll  = (h,h>^. 

Definition  2.1.  If  there  exists  a continuous  linear  functional 
G(u)  on  V such  that 


|E(u+h)  - £(u)  - G(u)h|  = o( ||h|l) 


(2.1) 


as  OhM  ♦ o,  then  the  linear  functional  G(u)  is  called  the  Frechet 
derivative  of  E at  u and  G(u)h  is  called  the  Frechet  differential  of 
E at  u with  increment  h. 

The  higher  derivatives  are  defined  in  a similar  manner.  Denote 
the  conjugate  space  of  V by  V*,  the  space  of  all  linear  functionals 
on  V.  Denote  the  norm  on  V*  by  ||*l|*. 

Definition  2.2.  If  there  exists  a continuous  linear  operator 
F(u)  from  V into  V*  such  that 


|lG(u+h)  - G(u)  - F(u)hl I • - o(llhll)  (2.2) 

as  llhll  ♦ 0,  then  the  operator  F(u)  is  called  the  second  Frechet 
derivative  of  the  functional  E,  and  E is  said  to  be  twice  differentiable. 
F(u)h  is  called  the  second  Fhechet  differential. 

Definition  2.3.  If  the  limit  ts  9^0  exists, 
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«(»,„)  . 11,  Si e^Ls  S<n> 
e ♦o  ® 


Ti.3) 

then  6E(u,h)  is  called  the  Qateaux  differential  of  E at  u with 
increment  h.  Similarly  for  the  second  Qateaux  differential  at  u with 
increments  h and  k, 

ffiCu.h.k)  = lim  SteVhP}.-  >E(“.h)  (2.4) 

X ♦O  A 

provided  the  limit  exists. 

From  the  above  definitions  it  can  be  seen  that  if  the  Frechet 

differential  exists  at  u,  then  the  Qateaux  differential  also  exists, 

and  the  two  differentials  are  equal.  Although  the  converse  may  not 

necessarily  be  true,  the  sufficient  conditions  are  provided  by  the 

following  theorem.  Die  proofs  of  the  following  four  theorems  may  be 

found  in  cooks  on  functional  analysis  such  as  references  [5],  [6]. 

Theorem  2.1.  If  6E(u,h)  exists  in  tlu-uoM  $11,  u > 0,  and  if  it 

is  uniformly  continuous  in  u and  continuous  in  h,  then  the  Frechet 

differential  exists  and  Q(u  )h  = 6E(u  ,h). 

o o 

From  the  viewpoint  of  studying  extremal  points  in  function  space, 
the  concept  of  Frechet  differential  is  essential,  but  from  the 
computational  standpoint  the  Frechet  derivatives  are  often  obtained 
through  Equations  (2.3)  and  (2.4)  whenever  the  conditions  stated  in 
Theorem  2.1  are  fulfilled. 

Since  Q(u)  is  a continuous  linear  functional,  by  the  Riesz 
representation  theorem  there  exists  an  element  V£(u)  -in  V*  such  that 


Q(u)h  = <VE(u),h> 


for  every  h I V.  VE(u)  is  oollod  tht  gradient  of  X at  u.  Because 
F(u)h  to  o continuous  linos r functional  on  V,  than 

Cr(u)hlh  - <Hj(u)h,h) 

where  Hj,(u)  is  a continuous  linear  operator  on  V.  Hg(u)  is  called  the 
Nessian  of  E at  u.  If  h is  a unit  vector,  <VE(u),h>  may  be  regarded 
as  a directional  derivative  of  E in  the  direction  of  h. 

Theorem  2.2,  Suppose  that  the  functional  E on  V has  a relative 
extremum  at  u*.  It  is  necessary  that  VE(u*)  ■ 0,  i.e.,  Q(u*)h  « 0 
for  all  h in  V. 

Theorem  2.3.  Suppose  that  the  functional  E on  V has  a relative 
extremum  at  u*  subject  to  constraints  <rfj(u)  « 0,  j«l,2,...,n.  Suppose 
that  V^(u)  exist  and  that  they  are  linearly  independent.  If 
VE(u*)  ft  0,  it  is  necessary  that  there  exist  unique  real  numbers 

^l’^2'***,^n  no*  *#ro  ®u®h  that 

n 

W„*>  (2.5) 

1 

Theorem  2.4.  If  the  Gateaux  differential  6E(u,h)  of  a functional 
E exists  at  each  point  of  some  convex  set  D C V,  then  for  any  u and 
u*h  in  D, 

E(u*h)  - E(u)  • &E(u+th,h)  (2,6) 

for  seme  t in  [o,ll,  and  similar  expression  holds  when  X is  an  operator . 

Definition  2.*».  The  operator  T mapping  V into  V is  called 
continuous  at  uq  t V,  if  for  any  sequence  which  converges  to  uq, 
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i.e.,  lim  ||u  -u  II  ■ 0,  the  itquraet F(u.)  ‘converges  to  F(u), 

n — * n ^ ° 

i.e.,  li®  ||r(u  ) - r(u  )|l  ■ o. 

n o 

n ♦ oo 


SECTION  III 

METHOD  OF  CONJUGATE  GRADIENTS 


3.1  Historical  Survey 

The  method  of  conjugate  gradients  was  originally  developed  for  solving 
linear  systems  of  algebraic  equations  independently  by  Hestenes  and  Stifel 
[7] , [8]  in  1952.  Hayes  [9]  extended  the  method  in  1954  to  solve  linear 
operator  problems  on  Hilbert  space.  Antosiewicz  and  Rheinbodlt  [lO|  in  1962 
gave  further  consideration  to  the  rate  of  convergence  of  the  method.  Fletcher 
and  Reeves  [ll]  in  1964  applied  the  conjugate  gradient  technique  in  minimizing 
positive  definite  quadratic  functionals  in  finite  dimensional  space.  Daniel 
[12]  in  1965  gave  an  improved  estimate  of  the  rate  of  convergence  and  discussed 
the  applicability  of  the  conjugate  gradient  method  to  nonlinear  operator 
equations.  In  the  area  of  application  of  this  technique  to  optimal  control, 
Lasdon,  Mitter  and  Warren  [l3]  , and  Sinnott  and  Luenberger  [14]  have  treated 
unconstrained  problems  with  considerable  success. 

3.2  Linear  Theory 

(a)  Sequence  of  Expanding  Subspaces 

Let  A be  a postive  def intie,  self-adjoint,  continuous  linear  opera- 
tor with  domain  V,  a real  Hilbert  space,  and  range  RC  v.  Then  there 
exists  a real  number  m such  that  (u,  Au)  m (u,u)  for  every  u in  V,  and 
A has  a continuous  inverse  A"1  whose  doatain  is  R and  range  V.  The  linear 
equation 


Au  - k 
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(3.1) 


has  a unique  solution  h ■ k~\  for  any  given  k t P*.  Suppose  u is  an 
estimate  of  h.  h-u  will  be  referred  to  as  the  error,  r ■ k-Au  will  be 
called  the  residual  of  u as  an  estimate  of  h,  and 

E(u)  * <h-u,A(h-u)>  (3,2) 

is  called  the  error  functional.  Since  A is  positive  definite  the 
problem  of  solving  Equation  (3.1)  may  be  treated  from  the  variational 
setting  by  minimizing  the  error  functional  in  (3.2). 

Finding  the  solution  to  Equation  (3.1)  by  minimizing  the  error 

functional  E in  Equation  (3*2)  using  iterative  procedures  often 

involves  a sequence  of  expanding  subspaces,  i.e.,  a sequence  of  closed 

linear  subspaces  B„  of  V such  that 
n 


mie  iterative  procedure  is  the  basis  of  the  conjugate  gradient 

method.  The  following  three  theorems  illuminate  the  underlying 

philosophy  of  the  Iterative  procedure. 

Theorem  3.1.  Let  B be  a linear  subspaee  of  V and  u in  V.  Then 

o 

the  functional  E in  (3.2)  aatisfies  E(uq)  £ E(uQ+y)  for  every  y in  B 
if  and  only  if  9E(uo)  or  the  residual  at  uq  is  orthogonal  to  B.  In 
particular,  E(uq)  £ E(u)  for  every  u in  V if  and  only  if  uQ  ■ h,  the 
solution  of  Equation  (3.2). 


The  assumption  that  A is  self-adjoint  is  not  essential  since  from  the 
theoretical  point  of  view  the  equations  Ax*k  and  Bx»b,  where  B»A*A 
and  baA*k,  A*  la  the  adjoint  of  A,  are  equivalent. 


Proof i Lot  y bo  any  non-sero  element  la  B.  Thon 
E(uo+y)  - *(u0)  ■ <y*Ay>  ♦ <vE(uo),y>  . 

Hoaoo  if  yl(u0)  la  orthogonal  to  B,  thon 
B(uo+y)  - I(uo)  - <y,Ay>  > 0 , 

or  E(u0)  1 I(u0+y)  for  every  y in  B.  If  for  every  y in  B,  thon  for 
•ay  rool  t 

E(uo+ty)  ■ E(\»o>  ♦ t<VE(u0),y>  ♦ t2<y,Ay> 

p 

which  implies  thot  t(yE(uQ),y>  t <y,Ay>  ^ 0.  This  oxprooolon  con  bo 

truo  for  oufflciontly  snail  t only  if  (VE(u  ),y>  ■ 0 or  yE(u  ) is 

o o 

orthogonal  to  B.  Finally,  if  u aiainiBoa  E on  V,  VE(u  ) « 2(Au  -k) 

o o o 

oust  bo  orthogonal  to  V,  and  honco  must  bo  sero.  But  A is  positive 

definite,  therefore  u « h. 

o 

It  is  interesting  to  consider  a geometric  interpretation  of  the 
statement  of  this  theorem.  E(u)  ■ constant  defines  a family  of 
ellipsoids  about  h,  and  the  gradient  of  E at  u is  orthogonal  to  the 
ellipsoid  through  u.  The  linear  subspace  B is  a hyperplane  through 
the  origin,  and  uq+B  is  a hyperplane  through  uq.  Suppose  it  intersects 
the  ellipsoid  E(u)  * E(uq)  ■ M.  Then  there  is  a region  on  the  hyper- 
plane within  the  ellipsoid  so  that  S(u)  £ N unless  uo*B  is  tangent  to 
the  ellipsoid  through  uQ  or  y£(uQ)  is  orthogonal  to  subspace  uq+B 
(see  Figure  3.1) . 

In  view  of  the  positive  definite  and  continuity  properties  of  E 
and  Theorem  3*1«  we  have  the  following  conclusion. 

Theorem  3«2.  Let  B be  a closed  aubspace  of  V.  There  exists  a 
unique  uQ  in  B that  minimises  E(u)  on  B,  and  (yE(uQ),y)  ■ 0 for  every 


Figure  3.1.  A geometric  interpretation  of  Theorem  3.1 

Theorem  3.3.  Let  ^ be  an  expanding  sequence  of  closed  sub- 
spaces  of  V,  and  v = U Bn*  L#tVun^be  the  of  points  such 

that  ur  c Br  and  E(un)  = min^E(u),  u t Bn*V.  un  «h  u n -*oo  . 

Proof:  Since  * BnJis  an  expanding  sequence,  ^E(un)  j is  a 
decreasing  sequence,  and  there  exists  a real  X so  that  E^)  * Xn  -»  x 
as  n -*oo.  E(un+ty)  *X  for  any  t and  y e V;  this  implies  that 

<y,A(h-un)>2  «£  (Xn-X)<y,Ay>  , 

so  that  for  y = u -u  . m/n 
n in 

[E(um)-Xl  ♦ [E(un)-Xl  i<Vum’A(Vu» 

But  A is  positive  definite,  which  implies  that^uR 'is  a Cauchy 

sequence.  Since  V is  complete,  therefore  u eh  as  n*oo. 

n 

Suppose  that  V is  separable  so  that  there  exists  at  least  one 
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linearly  independent  sequence  * Pn  in  V,  so  that  the  finite- 
diaensional,  and  hence  closed,  linear  subspaces  B spanned  by 

n Oo 

(p0.Pi Pn-1}  formed  an  expanding  sequence  with  U B ■ V.  As  < 

0 n 

practical  natter  in  applying  Theorem  3.2  we  must  solve  the  minimisa- 
tion problem 


for  each  n,  and  the  solution  will  be  expressed  as 


where  the  coefficients  or . depend  on  n.  It  would  be  convenient  for 


this  procedure  if  the  coefficients  would  be  independent  of  n,  and  this 
leads  to  the  following  topic. 

(b)  Conjugate  Direction  Method 


Definition  3.1.  let  p , p be  non-xero  elements  in  V.  If 


then  p^  and  pn  are  said  to  be  A-con jugate  or  A-orthogonal. 

The  iteration  method  in  which  the  sequence  of  non-xero  elements 
‘"Pg^that  is  chosen  to  satisfy  the  A-conjugate  condition  is  a 
conjugate  direction  method.  The  elements  ^p^p^,...^  may  be  determined 
before  tne  Iterative  process,  or  the  element  pQ  may  be  determined  at 
the  nth  Iteration.  Let  f t ) be  an  approximating  sequence  to 


* 


r 


\ ^0*^1' *** ,pn-l^ ’ <md  the  coefficients  a are 

IlJ 

independent  of  n. 

Leans  3.1.  Let^pn>be  s sequence  of  elenents  in  V whose  elements 
are  mutually  A-or thogonal . 1st  be  the  subspace  spanned  by 

{po,pl’***,pn-l")*  ^PP08*  th8t 

u « 0 
o 

Vl  “ VV. 

(-VKuJ.P,) 

then  minimizes  E(u)  on  B^. 

Proof ; It  follows  from  u^  « Vl^m-lVl’  ^*^^'pj> 

* (,E(Vi)'fj>*JVi<*Vi'',3)-  r°r  3 * "-1’  we  have  (VE(uB),Pj> 

* ^yE(um_l,pj^  ^ the  d#finltion  of  oB,  <7£(uj+1),pj>  . 0.  Hence 
(yE(un),Pj)  * 0 for  j*l,...,a-l  or  VEtu^)  is  orthogonal  to  Bb,  and 
the  assertion  follows  from  Theorem  3.1. 

It  is  interesting  to  observe  that  if  V is  finite  dimensional, 

say  n,  then  un*h.  *nd  the  iteration  always  converges  in  finitely 

many  steps.  Whenever  (VE(ub),pb>  * 0,  a uB  and  the  minimum  of 

E(u)  in  Bb  Is  also  the  minimum  in  Bb+^.  This  occurs,  for  example, 

when  u^ah.  The  assumption  that  the  Iteration  starts  with  uo*0  is  not 

essential.  For  if  u^O,  consider  the  problem  A(u+uq)  a k{  the 

iteration  uQaO,  u^^uh+[(-yE(um),pB>/2(pB,ApB>]pB  then  converges  to 
A 

the  solution  hah-uo.  As  an  immediate  consequence  of  the  above  lemma 
and  Theorem  3.3,  we  have  the  following. 


elements  are  mutually  A- orthogonal.  Let  Bq  ba  tha  eubspace  spanned  by 
fp  ,p. ,...,p  . } and  V ■ I Ib  . Lat  u ba  arbitrary,  and 


Then  u converges  to  the  solution  h 


(c)  The  Method  of  Conjugate  Gradients 


In  the  discussion  above  on  tha  conjugate  direction  method,  the 


determination  of  tha  sequence  of  vectors 


was  governed  only  by  the 


requirement  that  they  be  A-conjugate;  their  determination  remains 


relatively  arbitrary.  Frost  the  computational  standpoint,  it  is 


convenient  and  frequently  desirable  to  generate  the  p 's  at  each 


step  in  the  iteration  process.  The  method  now  introduced  is  the 


algorithm  used  by  Hestenes  that  generates  a particularly  useful  set 


of  conjugate  directions.  Each  direction  p is  generated  by 


"A-conjugate-iting"  the  gradient  vector,  and  thus  the  name  conjugate 


gradient  is  given 


The  iteration  is  defined  as  follows 


Having  obtained  uQ,  VE(un),  and  pR,  the  iteration  is  continued 
according  to  the  expressions  below. 


(3.3) 


<-VKu  ),p  > 
a « ■■  ” n 

” ^p-J- 

Vl  ■ "n  * V, 

5E<"„)  * «*Vi-W 

Wy.Ap,,) 

” <V*V 

Vi  ■ -'E(Vi>  * *fn 

Thtore.  lh.  quMtiti.e  d.fln.d  by  th«  iteration  proc.M 

above  satisfy  the  following  relations: 

^pm*Apn^  = 0 • ® / n 
(pa*vE(un)>  = 0 , a < n 

<pa’VE(un)>  s nVE(uB)|f2  , m£n 
<VE(un>,7E(un»  = 0 , a / n 

Proof:  From  the  defining  expressions  of  cr  and  g , the 

n n 

equalities  above  may  be  shown  by  induction. 

Corollary:  The  conjugate  gradient  aethod  is  a special  case  of 
the  conjugate  direction  method. 

Proof:  The  assertion  follows  from  expression  (3.4)  above. 

If  the  set  of  vectors  ^P0.P^»...^  generated  above  spans  V,  then 
the  solution  h would  be  achieved  by  Theorem  3.3.  However,  even  for 
the  cases  in  which  the  set  (po,Pr...}is  not  coe*lete  in  V, 
nonetheless,  we  can  still  claim  that  u^  converges  to  h.  This  is  an 
Important  and  desirable  fact  of  this  aethod.  The  next  theorem  is 


devoted  to  demonstrating  that  point. 

Lhm  3.2.  o_  as  dafinad  by  equation  (3*3)  is  boundad  abova  by 
UpBM2/2(pn,Apn). 

Proof:  Mp  M2  - |*E(u)|l2  ♦ I2  Jlpn  Jl2 

* n n n-x  n*i 


i I tVE(u  )|  I2 
n 


Tharafora, 


llpntl2  |tVE(un)M2  !lpnHi 

^Pn~Apn>  (P^.Ap^^  lh»icf„  i 


n*  pn  IMIOH' 

9 n 

IIP„U2 


” t l7S(u  )M* 

It 


-2o,n 


Lawn*  3.3.  Lat  a ■ sup^X.X  c a pact  rum  of  tha  poeltiva  dafinita, 
sal f-ad joint,  continuous  linear  operator  A Than 


I*  *n 


Proof:  Tha  leans  follows  from  tha  fact  that 


u.Au)  _ *pn,Apn* 


<pn’V 


Lawns  3.4.  X is  a strictly  decreasing  function  of  n,  i.a., 
X(uQ)  > I(un+1),  unless  tha  solution  ia  attained  at  tha  nth 
iteration. 


Proof: 


*S)*<Vl)  " I(un)-*(un+flrnpn) 

* ■2®n<pn»A(h‘'*n)>’<,S<pn'Apn> 
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p 


3 

(p  ,VE(u  )> 

«L 

* •j-  (pnt VE(un)) 


*. 

mT 


n 1 IVE(u  )[l2 


The  quantity  ■*£  ||yE(u  )H2  is  positive  unless  u =h,  which  is  the 
assertion. 

Theorem  3.6.  The  sequence  £un  ^obtained  from  the  above  conjugate 
gradient  method  converges  to  h,  the  solution  to  Equation  (3.1). 

Proof?  In  Theorem  3.3%  we  have  shown  that  the  sequence ^ufl^ 
converges  to  some  element  in  V,  say  u*.  0^un0 is  a monotonicaHy 

decreasing  sequence  that  is  bounded  below  by  zero,  hence  it  converges. 
It  follows  from  the  expression  in  Lemma  3.4 

lim  -y-  | |VE(u)||2  = lim  [E(u)-E(u .)] 
n -woo  2 n n -woo  n n+1 

= 0 . 

But  according  to  Lemma  3.3,  a Thus,  IIVE(u  ) I ! = | l2A(h-u  )|| 

no  n n 

converges  to  zero  as  n -w oo  . By  continuity,  A(h-u*)  = 0.  A,  being 
positive  definite,  implies  u*  * h. 


3.3  Nonlinear  Theory 

Aside  from  the  ease  of  linear  or  quadratic  functionals,  when  the 
Frechet  differential  of  the  functional  is  set  equal  to  zero,  the 
resulting  equation  is  nonlinear.  In  this  chapter  we  trill  consider 
the  extension  of  the  technique  discussed  previously  for  the  linear 


theory  to  solve  the  equation  of  the  font 


P(u)  » 0 

where  P is  a nonlinear  operator  napping  the  real  Hilbert  space  into 
itself.  Suppose  that  u*  is  the  solution  to  the  above  operator 
equation.  Then 

0 * P(u*) 

= P(u+(u*-u)) 

= P(u)+F(u)(u*-u)+o( I lu*-ul | ) 

where  F(u)  is  the  Frechet  derivative.  Therefore,  for  u sufficiently 
near  the  solution  u*,  i.e.,  Ilu*-ull  sufficiently  snail,  we  alnost 
have  a linear  operator  equation  with  the  linear  operator  F(u).  Since 
^F(u)  depends  on  u in  general,  if  we  generate  conjugate  directions  with 
respect  to  F(u),  we  can  at  nost  assert  that  any  two  consecutive  pn 
vectors  are  F(u)-conjugate,  while  the  other  vectors  are  approximately 
F(u)-con jugate  depending  on  how  near  u is  to  u*. 

Assume  that  the  error  functional  E (or  performance  functional,  as 
it  is  often  called  in  control  theory)  defined  on  V possesses  the  follow- 
ing representation  about  u 

E(u+h)  = E(u)+(P(u),h>  +|<F(u)h,h>  + OOlhlf3)  C3 

where  P(u)  and  F(u)  are  the  first  and  second  Frechet  derivatives, 
respectively.  Suppose  that  E attains  its  minimum  at  u;  then  by 
Theorem  2.2,  we  can  assert  that  the  linear  operator  F{u)  is  positive 
definite  in  seme  neighborhood  D about  u.  We  will  make  the  assumption 


that  P(u)  is  self-adjoint.  As  w«  ahall  aaa  in  tha  following  sections 
eoncaming  tha  application  of  tha  aathod  to  control  problems,  the 
operator  P(u)  la  indeed  linear,  positive  definite,  continuoua  and 
aelf-adjoint. 

In  view  of  the  results  above,  we  make  tha  extension  of  the 
algorithm  givan  previously  as  follows: 
let  uo  be  arbitrary 

G(u  ) » vE(u  ) 

O O 

po  = -G(uo) 

Having  obtained  un,  pn  and  VE(ua),  the  iteration  is  continued  accord- 
ing to  the  expressions  below: 


J = u + nr  p 
n+1  n n*n 


(3.6) 


where  ®n  is  the  smallest  positive  solution  a of  (G(u  +op  ),p  > = 0, 

n n rn 


G(Vi)  ■ *«Vi) 


a = 

n 


(3.7) 

(3.6) 


Vi  * -°(W  * paPn 


(3.9) 


To  determine  the  value  of  or  that  satisfies  the  equation 
(G(un+apn),pn>  s 0 is  a difficult  task  in  general,  but  we  will  make 
the  following  observations. 

Immma  3»5»  Let  D be  a convex  region  in  V containing  u*  such 
that  r(u)  la  positive  definite.  Then  E is  a convex  functional  or.  3. 
Proof:  Suppose  u^,u2  t D,  then  by  convexity  of  D, 


tu1+(l-t)u2  « D for  0 £ t £ 1.  The  rest  follows  immediately  from 
Equation  (3.5)  and  the  fact  that  the  weighted  sum  of  the  squares  is 
greater  than  or  equal  to  the  square  of  the  weighted  sums. 

Lemma  3.6.  If  E is  a convex  functional  on  a subset  W of  V,  then 
<^u:E(u)  £ M,  M being  a positive  constant^  C W is  convex. 

Proof:  E(tUj+(l-t)u2)  < tECu^+d+OECug) 

< M 

which  implies  the  assertion. 

As  the  consequence  of  the  above  lemmas,  we  have  the  following 
result. 

Theorem  3.7.  If  F(u)  is  positive  definite,  then  the  value  of  or 

that  minimizes  E(u  +op  ) coincides  with  the  value  of  or  that  satisfies 
n n 

<G(un+apn),pn>  = 0. 

We  have  thu6  reduced  the  problem  of  finding  the  solution  to 

(Q(u  +op  ),p  ) = 0 to  a one-dimensional  minimization  problem, 
n n n 

Theorem  3.8.  Suppose  that  F(u)  is  uniformly  bounded  and 
uniformly  positive  definite  in  Q,  the  closure  of^u:E(u)  £ E(uq)^  . 
For  the  sequence  ^ur^ generated  by  Equations  (3*6)  to  (3*9) « {g(uq)  ^ 
converges  to  zero. 

Proof:  Since  F(u)  is  uniformly  positive  definite  in  Q,  then 

there  exists  positive  constants  m and  M so  that  0 < ml  £ F(u)  £ MI, 

where  I is  the  identity  operator,  and  6 is  the  null  operator.  Let 

be  the  value  of  o that  minimizes  E(u  +crp _)*,  then  (G(u+orp),p)  * 0 

n n n n n 

and  consequently 


lft 


(Q(Wn)'Vl>  " (Q(Wn)-0(Wn)+Pn) 

- ||0(uB+ornpi|)|!2 

H (un+o*,n>^0  * -<Q(u„+0'nPn)’Pn> 

- |lQ(u  )||2 

n 

(«»♦<»„)  ■ <r<v«Vpn.p„> 

<Hllp ,|l* 
n 

whenever  un+apn  is  in  Q.  Since 

I!  * 0 

ON  Of 

n 

therefore 

*nM,,pnn2  > 

Firoo  Equation  (3»9)  *nd  the  fact  that  the  consecutive  elements  in  the 
sequence  r(u)-con jugate  thus 

«l|pnll2i  M||Q(ua)||2 

®n  ^ mA^  • iVom 

*(V<V  • K«D)«<a(u11).p11>  ♦ £ (o(un«Mi,n)Pii,pn> 

for  som  t t [0,1].  for  al  on 

*■  *<V-®I  •*  ♦ r or^Ml  ipnl  I2 

* *(V«I MV1 12  ♦?**!-  Ilatv11* 
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In  particular,  consider  or  ■ m/M,  then 


*‘Vl)  5 »„> 


* *(«„>-  Sj  llotun,|  l2  £j  1 10(»B)  1 1* 

i 1(^1  - ^ “»j  ||G(un)l!^ 

Therefore, ^J(un)^ is  ?.  monotonically  decreasing  sequence,  and 
hence  convergent  with  the  assumption  that  E(u)  ia  bounded  below. 
(«Vl>-£(»n>),°  as  n ■♦oo  implies  that  1 1 G ( ) 1 1 ♦ 0 as  n ♦oo  , 
As  s consequence  of  the  uniformly  positive  definite  property 
of  F(u)  in  Q,  and  supposing  that  Q is  compact,  we  then  have  the 
following  result. 

Theorem  3.9.  The  sequence  ^ur^ converges  to  a unique  element 
u*  in  Q,  the  solution  to  the  minimization  problem. 


3.^  Remarks 

(a)  In  finding  the  value  of  or  that  minimizes  E(un+apn),  we  may 
use  the  expression  or  given  in  the  linear  case,  namely, 


(5.10) 


as  a first  order  approximation  to  guide  the  initial  aearch.  The 
quantities  pft,  0(un>  and  F(un>  are  already  available  in  the 
computational  process,  thus  the  evaluation  of  or  does  not  involve  ouch 
work.  For  most  problems,  it  is  expected  that  this  approximation  would 
get  better  as  uQ  gets  closer  to  u*,  for  then  the  vectors  would  be 
closer  to  mutually  F(u*)-oon jugate. 
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(b)  Theorem  3.8  indicates  that  it  is  desirable  to  select  the 

initial  estimate  of  the  solution  u such  that  the  second  Frechet 

o 

derivative  of  E at  uq  is  uniformly  bounded  and  uniformly  positive 
definite.  We  may  use  this  fact  as  a guide  to  select  the  initial 
approximation. 

(c)  In  optimal  control  with  state  variable  constraint  applica- 
tions, the  set  of  admissible  controls  (this  term  will  be  made  precise 
later)  in  general  will  not  form  a linear  subspace  nor  even  meet  the 
convexity  hypothesis  in  the  discussion  above.  For  the  problem  in 
Section  VI,  the  admissible  control  set  possesses  the  necessary 
conditions  being  convex,  while  for  the  problem  in  Section  VII  it  does 
not.  The  modifications  made  to  obtain  convergence  in  computation  to 
the  desired  solution  are  presented  in  detail  there.  The  convergence 
is  along  the  expanding  sequence  of  sets  fl  Q^where  Q denotes  the 
set  of  admissible  controls. 

(d)  If  &n  is  set  equal  to  zero  in  each  step,  the  direction  of 

search  p^  would  be  along  the  negative  direction  of  the  gradient  of  E 

at  ur,  and  if  crn  is  selected  so  that  the  performance  functional  is 
minimized,  then  this  is  the  well-known  method  of  steepest  descent.  It 
is  worthy  to  note  that  in  the  method  of  steepest  descent,  the 
performance  functional  is  not  minimized  in  a sequence  of  expanding 
subspaces  as  it  is  in  the  conjugate  direction  methods. 

(e)  At  any  step  of  the  iteration  process,  we  can  start  anew  with 

only  a small  amount  of  labor  involved,  keeping  the  approximation  last 

obtained  as  the  initial  estimate. 

(f)  Other  variations  of  the  conjugate  gradient  algorithm  when 


is  Independent  of  u may  bo  found  in  tbo  poporo  of  Bootoaoa  and 
C73*  C81  yrhtn  tho  development  to  prooontod  la  groat  do  toil. 

SECTION  IV 

THE  CLASS  OF  CONTROL  PROBLEMS  TO  BE  CONSIDERED 

Our  ultimate  goal  is  to  apply  the  technique  develope  in  the 
previous  section  to  solve  the  class  of  control  problems  which  we  form- 
ulate below.  Supoose  that  the  dynamical  system  is  governed  by  the 
differential  equation 


*L. 


f(x,x) 


C4.1) 


where  x is  a real  n-vector  for  each  t,  called  the  state  of  the  system; 

u is  a real  m-vector  for  each  t,  called  the  control  vector;  and  f is  a 

real  n-vector  for  each  t,  that  is  twice  continuously  differentiable  in 

its  argwents  (t  will  be  interpreted  as  time  with  values  in  E14).  Let 

x(t  ) be  the  initial  state  of  the  system,  and  let  it  be  desired  to 
~ 0 

transfer  the  system  from  the  given  initial  state  to  some  final  state 
lying  on  a smooth  hypersurface 


♦ [x(tf)]  - 9 (4,2) 

where  the  terminal  time  t^  is  not  fixed,  while  the  states  are  confined  to 
whithin  a closed  region  in  En  given  by  the  inequalities 

gk(x)  k-l»2,,,,,N  (4.3) 

where  g^is  an  m-time  continuous  differentiable  function  of  x.  We  will 
callla  control  u an  element  in  the  Hilbert  space  of  piecewise  continuous 
functions  on  ft#  , tJ  with  inner  product  defined  as 


6 


B" 


denotes  an  m-dimenSional  Euclidean  space 
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fxt 

“ J u^Ou^Odt 


admissible  if  th#  corresponding  trajectory  in  E°  does  not  violate  the 
atate  constraints  above  for  all  t c Ct#,tf].  Denote  the  class  of 
admissible  controls  by  Q.  Let  the  performance  functional  be 


f f 

E(u)  « w[x(tf)]  ♦ L(u,x)dt 
J t ~ 


(4.4) 


or  alternatively  as  in  the  formulation  of  Mayer, 

E(u)  * *[x(tf)) 

a function  of  end  values  of  the  states,  where  x(tf)  is  an  augmented 
(n+1)- vector . In  the  following  x will  be  used  to  denote  either  the 
n- vector  or  the  augmented  (n+1 )- vector  without  further  specifying 
whenever  the  situation  is  clear  from  the  context. 

The  problem's  objective  is  to  find  the  control  u in  Q that 
minimites  the  performance  functional  while  satisfying  the  conditions 
(4.1),  (4.2)  and  (4.3).  Ve  will  make  an  assumption  that  there  exists 
a unique  solution  to  this  minimisation  problem. 


SECTION  V 


COMPUTATIONAL  CONSHBRATIONS 


$.1  Nomenclature 


For  each  control  u in  Q,  there  correaponda  a trajectory  in  e" . 
It  nay  conaiat  of  two  typea  of  area.  The  portion  of  a trajectory  in 


which  the  atatea  aatiafy 


will  be  called  an  interior  arc,  and  the  portion  that  aatiafies 


for  aone  k,  k*l,2,...,N  will  be  called  a boundary  arc.  A trajectory 


■ay  conpriae  entirely  a boundary  arc,  an  entirely  interior  arc,  or  a 
combination  of  interior  area  and  boundary  area  as  ahown  in  Figures  5.1 


Trajectory 


Constraining  surface 


Terminal  surface  t(x)  ■ 0 


Figure  5*1.  Trajectory  comprises  only  boundary  arc 


Flgur*  5.2.  Trajectory  comprises  only  interior  arc 


Figure  5«3>  Trajectory  comprises  interior  arcs  and  boundary  arc 
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The  trajectory  corresponding  to  the  control  u*  c Q that  minimises  the 
performance  functional  is  termed  the  optlaal  trajectory.  Consider 
tf  £ tQ,  i.e.,  tine  runs  forward.  The  smallest  value  of  t,  say  t,, 
for  which  x(t^)  lies  on  one  of  the  constraining  surfaces  g^(x)  * 0 is 
called  entering  tine,  and  xCt^),  the  entering  corner.  The  largest  t, 
say  t^t  for  which  xCt^)  lies  on  the  constraining  surface  g^Cx)  * 0 is 
called  exit  tine,  and  xCt^).  the  exit  corner.  For  simplicity  in  the 
discussion  below,  we  will  consider  only  the  cases  in  which  the  optimal 
trajectory  has  at  most  one  entering  corner. 

5.2  Control  on  the  Boundary  Arc 

For  the  period  t c [t^^]  along  the  boundary  arc,  the  states  are 
interrelated  by 

* 

gk(x)  « 0 (5.31 

It  follows  from  the  fact  that  along  the  boundary  arc,  the  constraint 
function  must  vanish  identically,  which  implies  that 


The  first  time  derivative  of  g has  a very  simple  geometric  interpreta- 
tion. It  states  that  the  boundary  arc  is  tangent  to  the  hypersurface 
g(x)  « 0,  or  normal  to  the  gradient  of  g.  That  is, 

If  • <**•  H>  rs.s) 

The  control  u will  be  determined  according  to  (5*5)  if  u appears 

^ If  two  or  more  constraints  are  involved  for  t t [t.,t_],  the  argument 
is  similar,  and  the  subscript  k on  g will  be  dropped  In  subsequent 
discussion. 
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I 

i 


explicitly  la  the  expression.  If  it  does  not,  we  asy  consider  the 
second  derivative  or  higher  derivatives  so  that  u will  appear  explicitly 

dN 

in  — 4f  « 0.  If  the  systen  is  controllable  [15],  the  existence  of  a 
dt" 

saallest  integer  N,  the  order  of  the  derivative  of  g for  which  u appears 
explicitly  is  assured.  From  (5.*0  in  particular  for  t « t^,  we  have 


g(x(t^)]  * 0 

(5.6) 

^ [x(t^)l  ■ 0 

j«l, ...,N-1 

(5.7) 

dtJ 


It  is  worthy  to  note  that  Equations  (5.6)  and  (5.?)  along  with  the 
control  u satisfying 


isply  that 


(5.8) 


(5.9) 


for  all  t i [tj.tg]. 

We  will  make  the  necessary  assunptlons  such  as  g has  no  singular 
point,  i.e.,  Vg(x)  i 0,  to  permit  a possible  unique  solution  for  u in 
terms  of  the  states  in  (5*8).  Actually,  we  need  only  to  sake  such 
assumptions  along  the  optimal  boundary  arc.  But  from  the  computational 
point  of  view,  in  particular  using  the  conjugate  gradient  technique, 
since  we  have  no  advance  knowledge  of  the  whereabouts  of  the  boundary 
arc  on  the  constraint  surf see,  the  above  provisions  are  necessary. 


5«3  lhe  Perturbation  Equations 

Suppose  that  u*  t Q is  the  optimal  control  of  our  minimisation 

problem,  and  that  x*  is  the  corresponding  optimal  trajectory.  Consider 
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u ■ u**8h  t Q,  wh ere  6 U * ml  numb er  and  h a piecewise  continuous 
function.  Lat 


x(t)  * x*(t)+Mt) 

ba  tha  trajectory  generated  by  u.  It  follows  from  (5.10)  that  wa 

dx(t)  dx*(t)  . d* 

SI TE  * • 


(3.10) 


(5.11) 


Since  x is  a solution  to  (4.1), 
&&  « f[x(t),u(t)] 


= f[x*(t)+9z(t),u*(t)+8h(t)] 

* f(x*(t),u*(t)l  a ||| ft  fe(t)  ♦ |||,0h(t)+o(e)  (S 

where  ||  is  tha  Jacobian  matrix  of  f with  respect  to  x,  ||  is  a row 
vector.  Tha  symbol  |,  indicates  that  the  quantity  is  to  ba  evaluated 
along  tha  optimal  trajectory.  It  follows  than  from  (5*10)  and  (3.12), 
and  ignoring  the  factor  o(8)  (this  will  not  alter  tha  ultimate  outcome 
when  the  limit  of  6 approaches  zero  is  taken),  that  z is  a solution  to 


(S.12) 


the  linear  differential  equation 


dz  dfi  _ . oi i . 


«(t0)  ■ o . 

On  the  boundary  arc,  the  control  and  states  are  further  subject  to 
N N 

H » 0.  Denote  by  Q. 
dt  it 


C5.13) 


■ ■iBmm 


oCx(t),u(t)3  ■ a[x*(t)+»*(t)tu*(t)+eh(t)] 


. G(x*,u*)+8  (|||„*>  * |^|,0h+o(0) 
Again,  Ignoring  the  factor  o(8)  we  have 


dQ  i 


* 0 


da  i 


If  f§L  * °»  then  it  follows  from  (5.13)  and  (5.15)  that 


dz  ar.  ori  rOG,  -.dGi  x 

TTF  3x***  “ Tu'*  W •»*> 

on  the  boundary  arc. 


df»  rdG^  .dGi 


(5.14) 


(5. IS) 


(S.16) 


5,i*  Firat  Frechet  Differential  of  the  Performance  Functional 
Assume  that  the  performance  functional  E as  defined  below 
satisfies  the  conditions  in  Theorem  2.1,  then  we  may  evaluate  its 
Prechet  derivative  by  formula  (2.3),  and  ultimately  obtain  the  gradient 
of  E.  Let  X be  a piecewise  differentiable  n-vector-valued  function  of 
t as  yet  unspecified.  We  will  call  X(t)  a costate  vector.  Let 

H(x,X,u)  a (X,f(x,u)>  (5.17) 

and  call  the  scalar  function  H the  Hamiltonian  of  the  system.  Treat- 
ing the  conditions  (4.1),  (4.2),  (5*6),  and  (5.7)  as  constraints, 
consider  the  performance  functional  E at  u a u*+8h  e Q. 

(*f 

E(u)  . ^Cx(tf)>#Cx(tf)Xu,S[x(t1)]>  + J [H(x,X,u)-(X,  |f)]dt  (5. 18) 

o 

where  u is  a constant  N-vector  and 
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sCxtt^l  » 


' gCxCtjH 


.N-l 


(5.19) 


for  convenience.  Making  expansions  about  the  optimal  trajectory, 

*[x(tf)]  - rf[x*(t*)]  ♦ «■  o(0)  (5.20) 

♦Cx(tf)]  * ♦Cx*(t*)l  ♦ ♦ o(0)  (5.21) 

S[x(t1)]  « S[x*(t*)]  ♦ ||l.t1dx(tl)  4 o(0)  (5,22) 

For  the  functional  I it  u#,  we  have 

fl 

E(u*)  « ^x*(t*)>t[x*(t*)V(u,S[x*(t*)>  + 1 [H(x*,X*,u*)-<X*,-g-)ldt 

Jto 

(*.23) 


Therefore, 


E(u)-E(u*)  = <|||..dx(t)>t^<|l|*,dx(t)>t^(l*,  |||.ti>dx(t1) 

♦ j f[H(x,X,u)-H(x*,X*,u*)-e<X,  «^>]dt+o(e)  (S.24) 

*o 

To  take  into  account  the  possible  discontinuities  at  the  entering  and 
exit  corners,  the  above  integral  will  be  written  as  three  integrals 
over  the  intervals  [ t0,t^),  (t^.tg)  and  (t^.tj.  ].  After  integration 
by  parts  is  performed,  (5*24)  becomes 
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E(u)-E(u*)  « <|||/dx(t)>t  ♦<||Kto(t)>t  +<u»  ||i#  dx(t1)> 

f • t •t1 

C1 

♦e<x(to),*(to)>-0(x(t^),r(t-)M)  C<|£,*>*<|§!  ,*>+|||h]dt 


♦e(x(tJ),*(tp>-ea(t2),z(t2)>+0 


ftf 

^0<X(tp,z(tp>-0(X(tf),z(tf)>+0  Jt+C<||,z>+<||Kz>+  |g|  hldt 


But  dxftj)  = 0z(tj)  + 0C6) 

where  i\  = tj-t*  , j=l,2,f 

and  z(t  ) = 0 

o 


(5.25) 


(5.26) 


dx(t.) 


■ <t£l  ♦&!  -'(tf).to(tf»*<X(tf),  -jj£-  V 


•t  *t 

zt 


dx(tp) 

♦<X(tJ>-X(tJ>,dx(t2»-<X(tJ>f  —3^  t2> 

dx(t“)  dx(t^)  jke  . 

+<x(t'),  — T2>-<x(tJ),  — • »i>*€L  H-X(t”)*X(t*» 


(fc-27) 


We  are  now  going  to  impose  the  following  conditions  on  the  eostate  X. 

(a)  Demand  that  X be  a solution  of  the  differential  equation 

4*  * * o*  cs.w 

dt  ox 

for  t c [t  ,t, ) and  t e (t_,t,]. 

0 1 c,  I 

(b)  On  the  boundary  arc,  for  t e (t^tj),  require  X to  satisfy 


dX  9f  . fit  .vrdG,-1 

dt  + 3SX  + <3u*x>W  ^ = 0 

(5.29) 

^ t^  and  tf  demand  that 

X(t‘)  = X(t *)  + ||l  » 

*1 

(5  50) 

<x’-§>  = <^i>t+ 

ii  ii 

(5.31) 

X(t")  *=  X(t+) 

(5.32) 

* <x,n>t+ 

t2  t2 

(5.33) 

x<  v * & v ♦ M<4f  > 

(5.34) 

J The  symbol  ( after  the  quantities  «r— , etc.  will  be  deleted 
henceforth.  • x 


Equations  (5.31)  and  (5.33)  indicata  that  the  Hamiltonian  is 
continuous  at  the  entering  and  exit  corners,  and  Equations  (5.32)  and 
(5*30)  show  that  the  costate  is  continuous  at  the  exit  corner  while  it 
is  discontinuous  at  the  entering  corner  with  a jump  equal  to  ^-1  u. 

°x  t 

With  the  above  conditions,  the  Frechet  differential  is  1 

ft 

G(u*)h  * ) ||  hdt  + J + ||  hdt  (5.35) 


and  consequently,  the  gradient  of  the  performance  functional 


VE(u)  = ||(x,X,u) 


(5.36) 


for  t e [to,t’)  and  t e (t*,^!. 


5*5  The  Second  Frechet  Differential  of  the  Performance  Functional 
Suppose  that  for  the  control  u » u*+0h  e Q,  the  corresponding 
trajectory  is 


x(t)  = x*(t)  + 0z(t ) + 02u)(t) 


(5/37) 


instead  of  (5.10),  where  z again  satisfies  the  differential  equation 
(5.13)  with  z(to)  = 0 and  (ii(tQ)  = 0.  In  view  of  (5.37),  let  us  now 
re-examine  the  expression  (5.24).  The  expression  (5.20)  becomes 

if! 

&X* 

(5.38) 

where 

, dx(t,) 

(5.39) 


2 

*[x(tfn  = ^x*(tp]+(||(tf),dx(tf)>+|<dx(tf),  2-|(tf)dx(tf»+o(02) 


? dx(tj 

dx(tf)  « 0z(tf)  ♦ 0 ®(tf)  + 0 ■ -jt-'  Tf 


j 1 


* 


Similarly,  expressions  (5.21)  and  (5.22)  beooae 


♦Cx(tf)l-tCx*(t*)V<|J(tf),dx(tf)>+  £<dx(tf),2-|<tf),dx(tf)>+o(02) 

bx 


(5.40) 


2 

S[x(t1)>S[x*(t*)>  ||(t1)dx(t1)+  |(dx(t1),^-|dx(t1)>+o(e2)  (5.41) 

The  integrals  in  (5*25)  become 

(*1  A ? dx(tT) 

- j (\,^(0*-»-0  u>))dt  s -0(X(t^),dx(t^)+  T^) 


\ ♦ 62  ^ ("j£,u))dt  ♦ o(02) 


(5.42) 


A P . . dx(tp) 

- j +(X,-j^(0»+02a))>dt  ■ 0<X(t2),dx(t2)+  — jg — t2> 

Jt2 

- 0(X(t^),dx(tg)+  T^*0  ^ +^H*z^dt'f®2  j ♦^3pu,^d* 

t.  t. 


♦ o(02) 


A P ♦ dx(tJ) 

■ ) +<X,«j^(0»+02a»)>dt»0(X(tJ),dx(t1)>+  ■ 


- 0(X(t‘),dx(t2)+ 


dx(t,) 


(5.43) 


t2>+0  ( +(^,s)dt+02  j 4,(j^,«u)dt'*,o(0i 
't,  t. 


(5.44) 


H(x,X,u)-H(x%X*,u*)-e<||,*>+e2(||,u»>+8(||,h> 


(5.45) 


From  Equations  (5.38)  to  (5.45),  in  view  of  the  conditions  imposed  on 
the  costate  X and  that  ^ = 0 along  the  optimal  trajectory,  we  have  for 
the  second  Frechet  differential 


(5.46) 

5.6  An  Approximation  for  the  Hessian 

gradient  of  the  performance  functional  follows  immediately 
from  the  first  Frechet  differential.  However,  some  further  steps  are 
necessary  from  (5.46)  to  obtain  the  Hessian  of  E which  is  to  be  used 
in  the  conjugate  gradient  algorithm.  Recall  that  * is  the  solution 
to  the  linear  time-varying  differential  equation 


for  the  unconstrained  region  given  by 
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r t 

s(t)  - L(t«tj)a(tj)  ♦ J L(  t , T ) |J(T)l^h(T)dT 

tj 


(S.«) 


where  L is  the  fundamental  matrix  and  t.  ■ t .tt.  s(t  ) ■ 0,  but 

j o c.  o 

z(t2)  as  yet  is  not  known.  Writing  ($.48)  as 


s ■ Th 


(5.48) 


and  letting 


■ra 


Formally,  we  have 


d2H 


&H\ 

'Sx'SZ 


(S.50) 


(5. SI) 


(*(tf),^^|  s(tf)>  * Th)t 

1 ix*  * 1 ax2  *f 


fVechet  differential  now  becomes 


(5.53) 


from  which  we  can  obtain  the  Heaaian  of  the  performance  functional. 

Since  (5.^7 ) describes  a time-varying  linear  system  and  the  initial 

state  is  either  at  rest,  *(tQ)  . 0,  or  *(t2>  which  is  small,  (when  u 

is  sufficiently  near  u*),  a "small”  h will  generate  a "small"  i (the 

converse  is  not  necessarily  true)  and  the  term  (h,^-|  h>  is  the  "most" 

5,2^,  Su4" 

dependent  on  h in  (5.53)  C 15l-  Furthermore,  — 5 is  positive  defini^ 
r - dud 

L15J«  and  therefore,  by  continuity  there  exists  a region  about  u*  for 

* 

which  — » is  positive  definite. 

Su* 

In  the  discussion  in  Section  III,  nieorem  3.8,  the  only  require- 

% 

ment  on  F(u)  was  that  it  must  be  uniformly  positive  definite  in  the 

neighborhood  of  u*  to  assure  convergence  of  the  iteration  process 

(assuming  also  that  the  controls  in  this  neighborhood  are  all 

admissible).  Hence  we  may  consider  constructing  a set  of  search 

directions  £pn  ^conjugate  with  respect  to  or  locally  conjugate 

A du 

to  be  precise  since  — j depends  on  u in  general.  This  simplification 
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provides  e considerable  reduction  in  copulation  tine  and  programing , 
and  «.  .ill  us.  Ms  approrination  belo.  in  our  copulations. 

5.7  Computation  for  the  Costate 

In  order  to  obtain  the  Hamiltonian  at  the  nth  iteration  on  which 

the  gradient  and  Hessian  are  based,  we  must  have  the  state  variable 

and  the  eostate  in  addition  to  the  control  u chosen  for  that  iteration. 

n 

Since  the  state  variable  is  continuous  and  the  initial  condition  is 
given,  solving  Equation  (4.1)  is  a straightforward  problem,  provided 
the  estimated  quantities  t^  t^  and  tf  are  settled.  We  will  elaborate 
this  point  in  the  next  subsection.  On  the  other  hand,  the  determina- 
tion of  the  costate  requires  more  considerations.  Since  the  boundary 
condition  for  the  costate  is  specified  at  the  terminal  time  tf  and  the 
costate  is  continuous  at  the  exit  corner,  thus  X(t)  for  t > (t*,tf] 
nay  be  determined  simply  by  solving  the  differential  equations  (5*28) 
and  (5*29)  backward  in  time  using  the  latest  estimated  control  and 
state  variable.  At  the  entering  corner,  when  t « t^,  the  costate  may 
be  discontinuous.  In  principle,  it  is  possible  to  determine  X(t“), 
it,  t^,  a total  of  N+n+1  unknown  quantities,  from  Equations  (5*30), 

(5*31) « (5*6)  and  (5.7)  ss  long  as  these  equations  are  independent. 

8ince  in  any  stage  of  the  Iteration  process,  the  time  at  which  the 
trajectory  reaches  the  constraint  surface  t^  is  in  general  not  equal 
to  t£,  hence  an  exact  solution  to  the  above  quantities  is  not  really 
essential  provided  that  some  means  are  taken  so  that  theae  quantities 
would  oonverge  to  the  desired  values  as  the  process  progresses. 
Initially,  a trial  and  error  technique  may  be  used  to  obtain  an 
approximation  to  these  quantities.  Depending  on  the  problem  at  hand, 
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frequently  iom  intuition  u to  tho  physical  nature  of  tho  problem  nay 
•orro  u a guide  to  the  guess  and  the  method  by  which  to  improve  the 
estimates  at  each  step.  Thia  la  the  Boat  difficult  part  of  the 
eoBputation  and  alao  one  of  the  aoat  time-consuming  portions  of  the 
iteration  process. 

After  the  estimate  of  l(t")  is  selected,  the  costate  in  the 
remaining  portion,  for  t c «*y  again  resort  to  solving  the 

differential  equation  (5.28). 


5*8  Entering  Time  and  Beit  Time 

Since  the  control  program  is  updated  at  each  step  according  to 


Vl  “ % * Vn 


(5.51*) 


the  new  trajectory  may  reach  the  constraint  surface  sooner  or  later 
than  the  previous  iteration.  In  other  words,  the  entering  tiae  in 
general  varies  with  each  iteration,  and  it  is  dictated  by  the  control 
chosen.  If  t^n^  is  larger  than  t^n**\  there  is  no  problem.  However, 
when  the  opposite  is  true,  then  some  extension  on  u^  oust  be  aade 
for  the  time  interval  (t^n\t^n+1^)  such  as 


W‘>  ■ 

or  aoae  convenient  extrapolation  based  on  u^t*0*)  and  the  rate  of 
change  of  un  near  t^n\  When  the  estimated  solution  is  near  the 
optiaua,  signify  by  relatively  small  values  of  VE(uR),  a sore  accurate 
determination  for  the  entering  tiae  and  the  entering  corner  being 


TTTaT 


denotes  tj  for  the  nth  iteration. 
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desirable  (also  the  terminal  time  and  terminal  state).  Some  refine- 
ment in  step  else  "dt"  in  solving  the  differential  equation  near  ^ is 
necessary  in  order  to  minimise  rounding  errors. 

Concerning  the  exit  time,  as  it  was  observed  by  McIntyre  and 
Paiewonsky  [l6],  the  conditions  for  leaving  the  constraint  surface 
cannot  be  translated  into  mathematical  statements  that  can  be  used  in 
a computing  process.  Again,  must  be  estimated  and  an  Improvement 
made  in  the  estimation  according  to  some  means  such  ns  to  increase  t2 
when  the  new  control  causes  the  trajectory  to  violate  the  constraint 
surface  and  decrease  t2  otherwise.  The  amount  of  suitable  change 
involved  depends  on  the  problem  at  hand.  Often  too  large  a change 
may  cause  the  trajectory  and  some  subsequent  trajectories  to  deviate 
greatly  from  the  optimal,  while  making  too  small  changes  would  waste 
unnecessary  computing  time. 


5.9  Determination  of  Optimum  or 


It  is  convenient  to  divide  the  state  variable  constraint  problem 
into  three  parts  in  the  following  discussions,  and  designate  them  as 
Region  I,  for  t c [t^t^?  Region  II,  for  t « and  Region  III, 

f or  t t (t2,tf].  For  Region  II,  the  computation  for  the  optimum  «n, 
the  value  of  or  that  minimises  the  performance  functional,  or  step  siae 
in  the  search,  is  not  involved  sinoe  the  control  on  the  boundary  are  is 
not  free  to  vary.  For  Region  III,  the  optimum  may  be  determined  by 
using  Equation  (3.10)  as  a guide  for  the  search  and  to  compute  the 
performance  functional  for  selected  values  of  or.  A quadratic 
interpolation  may  be  employed  to  improve  the  effectiveness  of  the 


••arch  and  to  reduce  tha  number  of  values  of  or  naadad  to  ba  considered 


The  computation  for  the  optima  a&  for  Region  I needs  further  attention 
first  of  all,  due  to  the  presence  of  the  constraint  conditions  (5.6) 
and  (5.7) » tha  values  of  or  to  ba  considered  oust  be  selected  in  such  a 


way  that  u + op  are  admissible  controls.  This  is  a one-dimensional 


minimization  problem  subject  to  some  side  conditions.  It  is  desirable 


to  limit  tha  number  of  values  of  or  to  ba  considered  so  that  the 


computational  time  is  reasonable  while  maintaining  a tolerable  accuracy 
on  the  approximating  solution  for  each  iteration.  Secondly,  the 


evaluation  of  the  performance  functional  is  not  as  simple  as  in  the 
case  for  Region  III  since  the  value  of  rfx(t,)]  will  not  be  known 


until  the  complete  trajectory  is  computed  which  includes  Region  III 


where  the  trajectory  is  as  yet  to  be  evaluated.  Some  equivalent 
condition  at  t < t.  instead  of  ^[x(t-)]  sometimes  may  be  used  as  in 


the  re-entry  vehicle  problem  below,  or  as  in  the  minimum  time  problem 
selecting  the  optimum  ern  to  be  the  one  that  is  nearest  to  the  value  of 
or  Equation  (5*10)  provided,  while  satisfying  the  condition  that 


u + or  p is  an  admissible  control 


The  following  diagram  shows  the  steps  in  the  computational  process 
Because  of  the  lack  of  advanced  knowledge  of  the  initial  values  for 
or  , I and  p , some  convenient  values  such  as  zeros  may  be  used. 


5.11  Substitution  of  Penalty  function  for  Constraints 


In  the  state  variable  constraint  control  problems  we  have 
discussed  above,  most  of  the  hardships  in  computation  arise  from  the 


constraint  requirements  (5.6)  and  (5*7)*  The  penalty  function 
teehnique  la  designed  to  alleviate  these  difficulties.  Instead  of 


attacking  the  problem  directly,  it  reformulates  the  control  problem 
with  state  variable  constraint  into  an  unconstrained  problem  wherein 


the  original  performance  functional  is  augmented  by  a non-negative 


penalty  term,  a function  of  the  state  variable  x which  increased  in 


value  with  trajectories  that  violate  the  state  variable  constraints 


Sty  selecting  a suitable  sequence  of  non-negative  penalty  functions  in 
the  iteration  process,  it  is  conceivable  that  in  many  eases  the 


desired  solution  for  the  original  problem  would  be  achieved  as  the 
limit  of  the  sequence  of  approximating  solutions  obtained  in  the 
iteration.  Indeed,  this  technique  has  been  given  rigorous  justifica- 
tions by  various  investigators,  Moser  [173,  Russell  [l8]  and 
Okamura  [19],  just  to  mention  a few.  For  most  penalty  functions  the 


intermediate  trajectories  usually  violate  the  constraints.  That  is 


some  portion  of  the  boundary  arc  is  approached  from  outside  of  the 


An  adaptation  of  this  technique  to  suit  the  conjugate  gradient 
computational  method  is  as  follows.  A new  performance  functional  is 


where  C is  the  trajectory  under  the  control  u and  the  non-negative 


penalty  function  n as  a function  of  x has  tha  properties  that  for  x 

within  the  constraint  set  tt  has  saall  values  relative  to  ^[x(tf )],  and 

for  x outside  of  the  constraint  set  n has  large  values  that  increase 

with  the  distance  (with  some  suitable  metric)  away  from  the  constraint 

surfaces.  And  as  a function  of  n for  a given  x,  n is  a monotonically 

increasing  function  for  x outside  of  the  constraint  set  and  conversely 

for  x within  the  constraint  set.  The  gradient  of  E*  in  general  does 

not  approach  zero  as  the  optimum  solution  is  near,  due  to  the  added 

penalty  term  I n(x,n)ds.  Therefore,  some  other  means  must  be  uqpd  to 
JC 

signal  that  the  optimum  solution  is  near  in  the  iteration  process. 
Comparison  of  the  values  for  ^[x(t^.)3  in  consecutive  iterations  often 
fail  whenever  the  performance  functional  has  a "flat  bottom"  feature. 
Often  direct  comparison  of  un+^  with  u^  is  necessary,  such  as  evaluat- 
ing the  quantity 

k 

**un+k“un^  * ^ an+j-l^pn+j-l*pn+j-l^  (5.56) 

J=1 

To  avoid  instability  in  computation  which  causes  the  intermediate 

♦ 

trajectories  to  swing  far  from  the  optimal  trajectory  and  may  sometimes 
cause  the  approximating  solutions  to  diverge,  the  penalty  function 
cannot  be  too  "harsh."  On  the  contrary,  the  solution  may  have  a very 
slow  convergence  rate  which  would  make  the  computation  inefficient. 

Some  compromise  must  be  made  so  that  each  iteration  brings  the 
approximating  solution  closer  and  closer  to  the  optimum  at  some 
reasonable  rate.  After  the  selection  of  the  penalty  function,  the 
computational  steps  are  the  same  as  the  one  given  above  in  Figure  5.4 


of  tho  blocks 


SECTION  VI 


A MINIMUM  DISTANCE  WITH  FORBIDDEN  REGION  PRQBLZH 


6.1  Problem  Description 


As  an  application  of  the  foregoing  discussion  to  ths  stats 


variable  constraint  problems,  let  us  now  consider  a problem  of  moderate 


computational  difficulty  so  that  the  features  of  the  conjugate  gradient 


method  can  be  observed  with  greater  clarity.  Suppose  that  among  the 
planar  curves  joining  the  point  (4,l/*0  and  some  point  on  the  parabola 
with  its  vertex  at  the  origin  while  avoiding  a circular  region  as 
shown  in  Figure  6.1,  it  is  desired  to  find  one  that  minimise  the  length 


of  the  curve.  The  control  version  of  this  problem  would  be  to  find 


the  time-optimal  control  for  a piecewise  smooth  path  satisfying  the 


specified  conditions  traversed  at  a constant  speed,  where  the  control 


Trajectory 


Terminal  surface 


Forbidden  region 


Figure  6.1.  Geometry  of  ths  minimum  distance  problem 


■ 


s 


variable  u is  taken  as  the  angle  foraed  by  the  tangent  to  the  path  and 
the  negative  x^-axia  (see  Figure  6.1). 

This  problem  will  be  aolved  using  the  conjugate  gradient  aethod 
two  ways.  First,  the  coaputation  will  be  carried  out  considering  the 
constraints  directly  and  then  eaploying  a penalty  function  to  convert 
the  constraint  problem  to  an  equivalent  unconstrained  one.  Finally, 
another  computational  technique,  the  popular  steepest  descent,  is 
studied  with  the  same  considerations  given  as  in  the  first  case  of  the 
conjugate  gradient  aethod. 

The  performance  functional  to  be  minimized  is 


E(u)  * /rf[x(tf)1 


‘ ■Xi{t!) 


- 


dt 


The  system  dynamics  can  be  written 

*1 


IT 


r -k  COS  U 


to2 

IT 


■ k sin  u 


(6.1) 


(6.2) 


where  the  constant  k will  be  taken  as  unity  in  the  sequel  for 
simplicity.  Letting  tQ  « 0,  the  initial  conditions  are 
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Ff 


xx(0)  « 4 

x2(0)  « 1/4  (6.3) 

x3(0)  * 0 


6.2  Formulation  for  Numerical  Computation  Using  Constraints  Directly 
The  Hamiltonian  associated  with  this  problem  is 

H(x,X,u)  = -X^  cos  u ♦ Xg  «in  u ♦ Xj  (6.4) 


and  therefore  the  costate  equation  along  the  Interior  arcs,  or  Regions 
I and  III,  is 

§ = O (6.5) 

2 

In  view  of  Equation  (5*34)  and  that  ^(x)  = -x^  and  t(x)  = x^+x|, 
the  terminal  time  tf 

x2(tf)  = 1 

X2(tf)  » 2x2(tf)  (6.6) 

x3(tf)  «=  -1 


According  to  Equations  (5.36)  and  (5*53) . we  have  for  the  gradient  and 
Hessian  of  the  performance  functional,  respectively, 

VE(u)  = X ^ sin  u + X2  cos  u (6.7) 

F(u)  a X^  cos  u - X2  sin  u (6.8) 


On  the  boundary  arc,  or  Region  II,  the  control  is  required  to 
maintain  the  trajectory  so  that  It  trill  lie  on  the  circle 

2 p 

(Xj-2)  + x2  ■ 1,  hence 


! j 
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u(t)  ■ CO**1  X2(t) 


or 


u(t)  ■ ■in“1[x1(t)-2]  (6.9) 

Ohing  Equation  ($.28),  tha  eoatata  aquations  for  Ragion  II  ara 

dX, 


« -cos  u(Xj^  ain  u ♦ X2  cos  u) 
d\ 

* ain  u(X^  ain  u ♦ X2  coa  u) 


(6.10) 


dX 

By  Equations  (5.J0)  and  (5«3l),  the  jumps  of  tha  eoatata  at  the  entering 
corner  ara  governed  by 

Xi<tp  ■ X^t*)  ♦ m 2[x1(t^)-2l 

X.(t")  ■ X-U.*)  ♦ u2x_(t^) 

21  21  21  (6.11) 
x3(t")  « x3(tj) 

dx.(t")  dx-(tr)  dx.(t*) 

W ”"3S  + W ~~R  * W ““JE 

dx2(t*) 

4 W — ir5- 

and  from  which 

X.(t!‘)[-xp(t.)+coa  u(0]  + X-(t*)[x. (t. )-2-sin  u(t7)] 

it  . -3_1 i-i 1 1— 1 -1,  — 1_  (6.12) 

2^Cxi(ti)-2]c°a  u(t1)-x2(t1)ain  u(t^)} 

It  is  worth  observing  that  in  tha  iteration  process,  precaution  must 
be  taken  to  avoid  overflows  in  computation  since  the  denominator  of 
(6.12)  may  vaniah  when  the  approximating  trajectory  is  tangent  to  the 
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r ' , 

circle  at  t^.  Whan  this  occurs,  tha  numerator  vanishes  also.  In  view 
of  tha  Uniting  processes  involved  in  Equation  (5*25),  we  nay  therefore 
apply  L'Hoapital's  rule  to  Equation  (6.12)  and  obtain 

U ■ [-X^(t^)sin  u(tj)  - X2(t*)cos  u(t^)] 

After  a control  is  chosen,  and  the  initial  conditions  (6.3)  are  given, 
the  differential  equation  (6.2)  can  be  solved  in  a straightforward 
manner.  To  be  able  to  evaluate  the  gradient  and  the  Hessian  of  E,  the 
costate  in  Regions  I and  III  is  needed.  The  values  of  the  costate  in 
Region  III  are  clear  from  (6.3)  and  (6.6).  By  solving  (6.10)  backwards 
from  tg  to  t^,  we  have  Xj Ct-^> , j=l,2,3«  Through  Equation  (6.11), 

Xj(t~),  j-1,2,3  nay  be  determined,  and  consequently  the  values  of  the 
costate  for  Region  I obtained. 

The  set  of  admissible  controls  Q consists  of  controls  that  produce 

trajectories  starting  at  (4,lA)  and  terminating  at  some  point  on  the 

2 p2 

parabola  x^+Xg  = 0 avoiding  the  region  (x^-2)  + x2  < 1.  If  v and  w 

are  admissible  controls  with  corresponding  trajectories  and  Cg  so 

that  the  forbidden  region  is  not  in  the  interior  of  the  region  D 

2 

bounded  by  curves  C^,  C2  and  x^Xg  * 0,  then  for  0 £ a £ 1,  the  control 
av+(l-a)w  would  be  admissible  since  the  resulting  trajectory  lies 
entirely  in  D.  Hence  Q is  convex,  and  according  to  Theorem  3*8,  if 
the  control  selected  for  the  initial  iteration  is  in  the  neighborhood 
of  u*  within  which  the  Hessian  of  E is  positive  definite,  the 
convergence  of  the  Iteration  process  is  assured. 

6.3  Formulation  for  Nunerical  Computation  Paine  Penalty  Function 
Let  us  now  axaalna  what  aodifieatlons  aunt  be  made  when  the 


10 


penalty  function  is  introduced  so  that  the  problem  with  state  variable 
constraint  becomes  an  unconstrained  problem.  Let 


n(x,n)  « .0l[(Xj-2)^x|lA^n^  (6.13) 

where 

f >+nf  1 i n < 10 

A(n)  « < l4+2(n-10),  10  £ n < 20 

i 34-t4(n-20),  n ^ 20 

As  n becomes  large,  the  contribution  to  the  performance  functional 
along  the  trajectory  exterior  to  the  circle  is  small,  and  n(x,n)  is 
positive  everywhere  except  for  one  point,  x^  « 2 and  « 0,  which  is 
zero.  Hence  it  possesses  the  desired  characteristics  stated  in  the 
previous  section. 

The  new  performance  functional  to  be  minimized  is 


E(u)  « -Xj(tf) 

ftf  f tf 

« - * dt  + I i?(x,n)dt  (6.1*0 

-»o  Jo 

The  equations  describing  the  system  dynamics  (6.2)  remain  the  same 
except  for  the  last  expression  which  becomes 

« 1 - n(x,n)  (6.15) 

The  initial  conditions  for  the  states  are  again  given  by  (6.3)* 
new  Hamiltonian  is 

H(x,X,u)  - cos  u ♦ \2  sin  u ♦ XjCl-«(x,n)l  (6*l6) 
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with  the  terminal  conditions  ss  given  by  (6.6).  The  gradient  and  the 
Heasian  are  again  given  by  (6.7)  and  (6.8).  Since  there  are  no 


constraint  conditions  involved,  the  computer  program  becomes 


considerably  simpler.  Any  piecewise  control  having  the  corresponding 


trajectory  satisfying  the  initial  and  terminal  conditions  is  an 


a&nissible  control,  and  Q is  again  convex 


6. A Computational  Results 


Many  different  estimated  controls  have  been  selected  for  the 


initial  iteration  in  the  course  of  studies  for  computation  as  a 


constraint  problem  and  as  an  unconstrained  problem  using  either  the 


conjugate  gradient  method  or  the  method  of  steepest  descent,  but  only 
one  is  presented  below  in  detail  for  each  way  of  solving  this  problem 
since  the  convergence  characteristics  are  similar  for  any  estimated 
control  for  which  the  iteration  process  converges. 


fronting  the  problem  as  a constraint  problem,  the  following 


initial  estimates  are  used 


• ■ 2.8  (see  figure  6.2) 
o 

The  change  in  9,  69  is  wdt  according  to  tho  following: 

(1)  If  tho  trajectory  under  tho  now  estimated  control  dooo  not 
violate  tho  constraining  clrolo  In  tho  kth  iteration,  thon 

-.1  , k < 8 

66  . ( -.01  , 8 £ k * 15 
-.0005,  k > 15 

(il)  If  tho  trajectory  undor  tho  now  ootinotod  control  violates 
tho  cona training  circle  as  shown  below,  whore  A is  tho  point  where  the 
trajectory  loaves  tho  circle  in  tho  (k-l)st  iteration  and  B is  tho  last 
intersection  of  tho  trajectory  undor  tho  now  control  and  the  circle, 
thon 

69  ■ length  of  arc  AB 


* W 


Treating  the  problem  as  an  unconstrained  problem,  the  initial  estimate 
was 

u (t)  « 1 - .4t 
o 

The  fourth  order  Runge-Kutta  method  was  used  in  solving  the 
differential  equations  with  step-size  dt  » .01,  except  for  the  small 
intervals  slightly  before  the  trajectory  intersects  either  the  terminal 
curve  x^+x^  * 0 or  the  constraining  circle,  where  dt  «=  .0002. 

The  computed  results  are  shown  in  Figures  6.3  to  6.11. 

Comparisons  are  made  between  the  different  approaches  of  solution 
whenever  possible.  In  the  following,  whenever  there  is  no  mention  as 
to  whether  the  solution  is  obtained  by  using  constraints  directly  or  by 
using  a penalty  function  it  is  understood  that  the  first  is  used. 
Figures  6.3  and  6.4,  respectively,  show  the  approximating  controls 
uR(t)  for  various  n computed  by  the  conjugate  gradient  method  and  the 
corresponding  trajectories.  The  values  of  the  performance  functional 
vs.  n are  given  in  Figure  6.5,  and  the  gradients  of  the  performance 
functional  as  a function  of  time  for  various  n are  shown  in  Figures 
6.6  and  6.7.  Figures  6.8  and  6.9  provide  some  comparisons  of  the 
conjugate  gradient  and  steepest  descent  methods.  Since  the 
convergence  characteristics  of  the  approximating  controls  and  the 
corresponding  trajectories  in  Region  III  depend  mainly  on  the  choice 
of  ®n  or  the  exit  corner,  only  their  convergence  characteristics  in 
Region  I are  considered.  Figures  6.10  and  6.11,  respectively,  show 
the  approximating  controls  uR(t)  for  various  n computed  using  the 
penalty  function  approach  and  the  corresponding  trajectories. 


Approximating  Controls  Conputed  by  Method  of  Conjugate  Gradients 


Gradient  of  performance  functional 


flguro  6.9.  Deviation  of  approxiaatiac  control* 
fro*  optlaal  by  method  of  conjugate 
gradlanta  and  aethod  of  ataapaat 
doaeont 


oxiaating  controls  computed  by  Method  of  conjugate 
gradients  using  penalty  function 


Ioj*uo3 


6.11.  Trajectories  computed  by  method  of  conjugate  gradients 
using  penalty  functions 


SECTION  VII 


MAXIMUM  RANGE  OF  A REENTRY  VEHICLE  WITH  SKIP  ALTITUDE  CONSTRAINT 
7.1  Problem  Description 

Suppose  it  Is  desired  to  obtain  an  optimal  angle  of  attack  program 
for  a nonthrusting  reentry  vehicle  so  that  its  range  is  maximised  for 
some  given  initial  velocity  and  position  while  not  exceeding  a certain 
altitude  limit  once  it  has  gotten  below  the  limited  altitude.  The 
altitude  limitation  occurs  in  many  practical  situations  whenever 
maneuverability  of  the  vehicle  is  required  since  at  high  altitudes  the 
air  is  too  thin  to  provide  sufficient  aerodynamic  forces.  The 
parameters  used  in  this  problem  will  be  selected  to  resemble  a typical 


real  situation. 


Th*  earth  will  be  considered  as  spherical  and  non-rotating.  Thj« 
state  variable  constraint  control  problem  may  be  formulated  as  follows 
Let  Xj  be  the  range  along  the  earth;  »2  be  the  altitude  of  the  vehicle 
Xj  be  its  velocity;  and  be  the  flight  path  angle  relative  to  local 
horizontal.  The  constraining  surface  is 


The  performance  functional  to  be  minimized  is 


where  tf  is  the  time  when  the  vehicle  is  at  an  altitude  of  70,000  feet 
The  motion  of  the  vehicle  is  given  by  the  following  equations: 


P.  the  atmospheric  tensity  « .00238  sxp  (-x^^.ooo)  slugs/ft5 
R,  the  radius  of  aarth  * 2.1  x 107  fast 
S * tha  effective  surface  in  ft2 
m » tha  mass  of  tha  vahicla  in  slugs 


Let  tQ«0  and  the  initial  conditions  be 
3^(0)  = 0 feet 
x2(0)  = 340,000  feet 
Xj(0)  = 28,000  feet/seeond 
X|f(0)  = -.14  radians. 

7.2  Analysis  for  Numerical  Computation 

The  Hamiltonian  associated  with  this  problem  is 


(7.5) 


H(x,X,u)  = X1 


3 


cos  + XgXj  sin  x^ 


- Xj[(.2?4  + 1.8  sin2u)px|  ^ + g sin  x^] 


+ X,, 


.6  sin  2upx3  ^ ^ cos  x4  + cos  x^ 


(7.6) 


and  therefore  the  costate  equations  along  the  interior  arcs,  or  Regions 
1 and  III,  are 
dX 

-fir-  = Xxx?  cos  x4 — 5 - X3  [(.274  + 1.8  sin2u)xf  p/24,000 

(R+Xj) 

+ 2g  sin  xlf/(R+x2)]  + X^.6  sin  2ux^  p/24,000 

(7.7) 


2g  cos  x^  x,  cos  x^ 

xp*p  - 


(R+x2)* 


66 


In  view  of  Equation  (5.34)  and  that  ^(x) 


at  the  terminal  time  t 


According  to  Equations  (5*36)  and  (5»53)  we  have  for  the  gradient  and 
Hessian  of  the  performance  functional,  respectively, 


where  an  approximation  has  been  made  for  the  Hessian  as  discussed  in 


Section  V 


On  the  boundary  arc,  or  Region  II,  according  to  (5*4),  we  have 


Equation  (7.12)  implies  that  Xj 
Equation  (7.13). 


0,  and  consequently,  it  follows  from 


Therefore,  the  control  along  the  boundary  arc  is  given  by 


and  the  motion  is  given  by 


It  follows  from  Equation  (5.28)  that  the  costate  equations  for 
Region  II  are 


•"«■**** 


0 


dxl 

IT 


dX2 

IT 


Xlx3 


R 


(R*x. 


■p  - XjC(.274  ♦ 1.8  ain2u)x2^j  o/24.000] 


2ux?  o/24 , 000  - 

♦ (Xj  1.5  tan  2u  - X^/x^)[x2 

♦ x2/(R^x2)2  - 2g/(R+x2)] 


- X, 


.6  sin 


0 *6  sin  2u/24,000 


(7.17) 


* -X^/d-^Xg/R)  ♦ Xj(.274  ♦ 1.8  8in2u)ox^ 

- X^C.6  sin  2up  ♦ g/x2  ♦.  l/(R«-x2)] 

- (Xj  1.5  tan  2u  - X^/x^ 


) 


) 1.2  sin  2upx 


3 "5b 


dX^ 

IT  s "X2x3  * X38 

By  Equations  (5.30)  and  (5.31)  *nd  tha  fact  that  x^(t^)  « 0,  the  jumps 
of  the  costate  at  the  entering  corner  are  governed  by 
xi(t“)  * x1(tp 
X^tt^)  * uq  ♦ X2(tJ) 

x3(t")  « X3(tJ)  (7.18) 

X^(t“)  « i^Xj  ♦ X^(t*) 

dx,  dx.  dx,  dx, 

X3  It  1 1 - * X4  *H.  - * X3  It  • ♦ * X4  It"1  + 

1 *1  *1  *1 
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It  follows  fro*  ths  last  two  expressions  that 


[dx,  dx_ 

w it  1,4.  - . 

L h t 


Si 


1J 


yv  irl  - 


v*;> 

^7 


(7.19) 


u as  a function  of  the  state  and  costate  vectors  is  not  explicit,  and 
o 

a trial-and-error  procedure  such  as  the  one  described  in  Section  V will 
be  used,  From  the  pnysics  of  the  problem,  it  is  clear  that  in  order  to 
minimize  the  drag  force  (the  first  term  on  the  right  side  of  Equation 
(7.*0),  and  consequently  the  greater  range  the  vehicle  would  travel, 
it  is  desirable  for  the  vehicle  to  make  its  flight  at  the  limit 
altitude  oi  200,000  feet  for  some  period  of  time,  provided  that  it  has 
sufficient  energy  to  return  to  that  altitude  after  it  has  been  below 
200,000  feet  once.  This  fact  will  serve  as  a guide  in  selecting  the 
initial  estimate  of  the  control  in  Region  I. 

The  continuity  property  of  the  state  variable  implies  that 
x4(t-)  = 0 which  in  turn  implies  that  the  admissible  controls  must 
have  their  corresponding  trajectories  tangent  to  the  constraining 
surface  at  the  entering  corner.  Consequently,  the  set  of  admissible 
controls  Q does  not  possess  the  same  property  as  the  problem 
considered  in  Section  VI,  where  the  admissible  set  of  controls  is 
convex.  Moreover,  it  is  true  that  every  neighborhood  of  an  admissible 
control  has  at  least  one  non-admissible  control.  We  can  justify  this 
statement  heuristically  as  follows:  If  u is  an  admissible  control, 
then  we  can  select  a control  v so  that  v(t)  » u(t)  for  t c [o,t^-6], 
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6 > 0,  and  for  t > t^-6  select  ▼ bo  that  th#  corresponding  trajactory 
is  not  tangent  to  tha  constraining  aurfaca  (i.a.,  ▼ la  not  an 
admissible  control).  Sinca  6 is  arbitrary,  ||u-v||  say  ba  made 
arbitrarily  snail. 

From  tha  computational  viewpoint  it  is  sonatinas  advantageous , 
such  as  in  this  problem,  to  relax  the  requirement  that  the  estimated 


control 


for  each  iteration  be  admissible.  In  other  words 


eonputational  tine  nay  be  considerably  reduced  if  the  optimal  control 


u*  is  approached  along  a "path"  whose  intermediate  "points"  may  be 


nonadmissible  conceptually  as  shown  in  Figure  7.2 


The  set  of 

admissible  controls  Q 


Figure  7.2.  Illustration  of  the  path  leading  to  optimal  control  u 


Let  us  now  make  this  statement  specific.  Let  the  modified  entering 


time  t.,  be  the  smallest  t that  satisfies  either  one  of  the  following 


two  conditions 


The  flight  path  angle  vaniahea  also  before  the  pull  up 


The  adaptation  of  these  two  conditions  to  the  computer  program  is 
simple.  Denote  the  set  of  admissible  controls  for  the  kth  iteration  by 
Q^,  and  let 

\ = {ujx^)  < and  200,000  - x^)  < 5,00Q/k}  (7.20) 

The  sequence  of  admissible  control  sets  as  defined  above  approaches  Q 
as  k*»  oo, 

It  is  worthy  to  note  that  in  Region  I the  problem  of  determining 
or  that  minimizes  the  performance  functional  is  equivalent  to  the 
problem  of  maximizing  x^  for  a fixed  Xj^  in  Region  II,  thereby 
eliminating  the  computation  of  the  trajectory  in  Region  III  on  which 
the  performance  functional  is  based. 

7.3  Computational  Results 

The  initial  estimate  of  the  controls  for  Regions  I and  III  and 
the  exit  time  are 

UQ(t)  for  t e [0,^1  as  shown  in  Figure  7.4 
uQ(t)  = .42  radians,  t e [t2,tf] 
t^0^  = 235  seconds. 

The  change  in  t2,  At  is  made  according  to  the  following: 

(i)  If  the  trajectory  under  the  new  estimated  control  does  not 
violate  the  constraining  surface  x 2 - 200,000  = 0 in  the  kth  iteration, 
then 

r-2  k < 4 

At  » / -1  5 £ k S 8 

V.-.3  k>8 
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(ii)  If  the  trajectory  under  the  new  estimated  control  violates 
the  constraining  surface  as  shown  below  where  A is  the  exit  corner  in 
the  (k-l)st  iteration  and  B is  the  last  intersection  of  the  trajectory 
Ck  under  the  new  control  and  the  constraining  surface,  then 


dlatane*  between  A and  B 


Figure  7.3.  trajectory  violating  constraining  surface 

The  fourth  order  Runge-Kutta  method  was  used  in  solving  the  state 
and  costate  differential  equations  with  step  sizes  dt  = 1 for  Regions 
II  and  III  and  dt  * .5  for  Region  I,  where  greater  accuracy  was 
demanded,  except  for  the  small  intervals  slightly  before  the 
trajectory  intersects  either  the  terminal  surface  or  the  constraining 
surface  where  dt  * .01. 

Figures  7.4  and  7*5,  respectively,  show  the  approximating 
controls  uQ(t)  for  various  n computed  by  the  conjugate  gradient 
method  for  Region  I and  Region  III.  As  expected  from  physical 
consideration,  the  control  program  for  the  vehicle  when  it  is  above 


W 
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density  is  extremely  small.  During  the  ten  iterations,  no  significant 


changes  occurred  in  u^(t)  for  the  first  30  seconds,  which  corresponds 
to  x^  greater  than  220,000  feet,  consequently  the  curve  representing 

for  t < 30  may  contain  considerable  amount  of  uncertainty.  This 
fact  may  be  seen  by  examining  the  Hessian  of  the  performance  functional 
in  Figure  7.7  which  has  relatively  small  values  for  t < 30,  since  the 
Hessian  may  be  interpreted  as  weighting  factors  in  the  equations  for 
the  determination  of  or  and  8 in  (3. 1C)  and  (3.8).  The  Hessian  ha6 
relatively  large  values  in  the  vicinity  of  t = 55  which  indicates  that 
the  control  program  on  the  optimal  trajectory  during  that  period  of 
time  is  critical.  The  performance  functional  vs.  n is  given  in 
Figure  7*6,  and  it  indicates  that  after  the  8th  iteration  the  changes 
in  the  performance  functional  become  very  small  in  comparison  with  the 
changes  occurred  in  the  earlier  iterations.  The  computed  results  for 
all  the  states  after  ten  iterations  are  presented  in  Tables  7.1,  7.2, 
and  7-3,  and  the  costates  are  shown  in  Figure  7.8.  The  computed 
results  show  a very  close  agreement  with  the  optimality  relationship 
that 

u*(t)  tan"1  (3^) 

1 I 

as  required  by  the  Weierstrass  condition. 

| 

/ 

/ 





74 


night 


Time 

(seconds) 

Range 

(miles) 

Altitude 
(105  feet) 

Velocity 
(105  ft/sec) 

Angle 

(degrees) 

! 0 

0.0 

340.00 

28.000 

-8.0214 

5 

22.458 

320.53 

28.020 

-7.9626 

10 

44.954 

301.19 

28.037 

-7.9006 

15 

67.487 

281.99 

28.050 

-7.8319 

20 

90.049 

262.98 

28.051 

-7.7491 

25 

112.63 

244.20 

28.031 

-7.6365 

50 

135.21 

225.79 

27.967 

-7.4622 

35 

157.73 

207.99 

27.818 

-7.1647 

40 

180.12 

191.30 

27.520 

-6.6389 

45 

202.23 

176.50 

26.994 

-5.7397 

50 

223.85 

164.69 

26.191 

-4.3533 

55 

244.78 

156.89 

25.175 

-2.5466 

60 

264.89 

153* ^7 

24.127 

- .64277 

65 

284.21 

153.86 

23.227 

.96062 

70 

302.86 

156.95 

22.544 

2.0435 

75 

321.04 

161.54 

22.053 

2.6041 

8o 

338.86 

166.70 

21.699 

2.7499 

85 

356.43 

171.78 

21.428 

2.6317 

90 

373.79 

176.51 

21.212 

2.4488 

95 

391.00 

180.83 

21.034 

2.2374 

100 

408.07 

184.72 

20.885 

2.0100 

105 

425.03 

188.15 

20.757 

1.7740 

no 

441.90 

191.14 

20.647 

1.5338 

115 

458.68 

193.68 

20.548 

1.2920 

120 

475.38 

195.78 

20.460 

1.0500 

125 

492.01 

197.43 

20.379 

.80862 

130 

508.58 

198.65 

20.304 

.56854 

135 

525.09 

199.45 

20.233 

.33703 

140 

541.55 

199.86 

20.167 

.13610 

144.4  555*98  199.96  20.111  .00001 


TABLE  2 

The  dynamics  of  the  maximum  range  problem 
after  10  iterations,  region  II 


C(t-t^)  seconds]* 


0 

5 

10 

15 

20 

25 


50 

35 

40 

45 

50 

55 

60 

65 

70 

75 

79.47 


5 

6 
8 


700.55 

716.25 

731.88 

747.44 

762.92 

778.33 

793.66 

806.40 


20.111 

20.025 

19.938 

19.852 

19.765 

19.677 

19.590 

19.501 

19.413 

19.323 

19.233 

19.143 

19.051 

18.959 

18.866 

18.773 

18.693 


Hie  entering  time  t^  = 144.4. 
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TABLE  3 


(seconds) 


628.75 


The  dynamcis  of  the  maximum  range  problem 
after  10  iterations,  region  III 


0.0 

148.50 

289.55 

422.13 
544.26 
654.91 
754.10 
841.25 

915.13 
975.31 

1022.1 

1055.9 

1077.6 

1085.1 


Altitude 


200.00 


198.90 

192.94 

183.95 
177.35 
172.87 
165.86 
156.10 
146.37 
135.78 
121.74 
103.47 

82.420 

70.000 


18.693 

17.776 

16.828 

15.663 

14.277 

12.865 

11.446 

9.8873 

8.2082 

6.5411 

4.9259 

3.3970 

2.0024 

1.3706 


ight  Path 


(degrees) 


0.0 

- 0.20551 

- 0.56930 

- 0.62105 

- 0.38409 

- 0.45921 

- 0.88688 
- 1.1610 

- 1.3492 

- 2.1070 

- 3-7584 

- 6.8382 
-12.231 
-19.150 


Exit  time  t2  « 79.47  + 144.4 


Range  at  exit  time  R_  = 806.40  miles 


SECTION  VIII 


CONCLUSIONS  AND  REMARKS 


From  the  computational  experience  baaed  on  the  problema  studied 
above,  the  method  of  conjugate  gradients  has  been  shown  to  be  a useful 
computational  tool  in  solving  both  linear  and  nonlinear  optimal  control 
problems  with  state  variable  constraint.  The  method  is  basically 
simple  and  relatively  easy  to  program.  Although  the  search  directions 
are  only  locally  conjugate  with  respect  to  the  second  Frechet 
derivative  of  the  performance  functional,  they  still  provide 
satisfactory  convergence.  The  results  presented  in  Section  VI  indicate 
that  the  conjugate  gradient  method  has  a higher  rate  of  convergence 
in  comparison  with  the  method  of  steepest  descent,  but  the  difference 
in  the  rate  of  convergence  is  less  pronounced  for  this  constraint 
problem,  as  compared  with  the  cases  of  unconstrained  problems  reported 
by  other  investigators  [l3l,  [l*0,  because  of  the  following  reasons: 

(i)  The  set  of  admissible  controls  Q is  restricted,  and  consequently 
only  small  step  size  in  the  search  direction  is  permitted  in  Region  I. 
That  is,  the  convergence  is  along  the  expanding  sequence  of  sets 
< b n <0  instead  of  expanding  sequence  of  subspaces,  (ii)  The  rate 
of  convergence  in  Region  III  depends  heavily  on  the  choice  of  the  exit 
corner  in  each  iteration.  A considerable  portion  of  computational 
time  in  each  iteration  is  devoted  to  the  determination  of  the  optimum 
step  size  in  the  search  (although  the  exact  optimum  is  not  essential) 
and  the  determination  of  jumps  in  the  costate  at  the  entering  corner. 


To  assure  that  the  sequence  of  approximating  controls  converges 
to  optimum,  it  suffices  to  have  the  initial  estimated  control  so  that 


I 


the  ••cond  Frechet  dexivative  evaluated  there  is  positive  definite. 
The  sethod  of  conjugate  gradients,  like  many  other  optimisation 
techniques,  cannot  differentiate  local  minimum  from  absolute  minimum, 
and  consequently  the  initial  estimated  control  must  be  selected 
cautiously  unless  the  given  problem  is  known  to  have  only  one  minimum . 
In  converting  the  constraint  control  problem  to  an  equivalent 
unconstrained  one  by  introducing  a penalty  function,  the  computational 
process  involves  more  time  in  contrast  to  the  approach  which  considers 
the  constraints  directly,  but  it  requires  less  programming  work.  Its 
effectiveness  depends  heavily  on  the  proper  choice  of  the  function  n. 
We  have  treated  only  the  control  systems  that  are  time-invarying,  but 
the  extension  of  the  conjugate  gradient  method  to  encompass  the  time- 
varying  systems  is  straightforward. 
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